/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package Test;

/**
 *
 * @author Gebruiker
 */
import Jama.EigenvalueDecomposition;
import Jama.Matrix;
import java.util.Random;

/**
 * This class represents the blind source seperation procedure. 
 * @author Cliff
 */
public class BSSTest {
    public BSSTest(double [][]source){
        Matrix bBarMatrix = new Matrix(source.length, source.length);
        Random random = new Random();
        for (int i = 0; i < bBarMatrix.getRowDimension(); i++) {
            for (int j = 0; j < bBarMatrix.getColumnDimension(); j++) {
                bBarMatrix.set(i, j, random.nextDouble());
            }
        }
        this.random = bBarMatrix;
    }
    private Matrix whiteningMatrix;
    private Matrix random;
    /**
     * Performs prewhitening, necesarry for performing BSS.
     * @param source One signal with all channels.
     * @return 
     */
    public double[][] preWhiten(double[][] source) {
        int m = source.length, channels = source.length; //channels
        int n = source[0].length, samples = source[0].length; // samples
        Matrix muVector = new Matrix(channels, 1);
        Matrix covarianceMatrix = new Matrix(source.length, source.length);

        //Calculate mu vector
        for (int chan = 0; chan < channels; chan++) {
            double sum = 0;
            for (int sample = 0; sample < samples; sample++) {
                sum += source[chan][sample];
            }
            double avg = sum / samples;
            muVector.set(chan, 0, avg);
        }

        //Center the data. From every element, the mean of that channel is substracted.
        double[][] result = new double[m][n];
        for (int i = 0; i < m; i++) {
            double substraction = muVector.get(i, 0);
            for (int j = 0; j < n; j++) {
                result[i][j] = source[i][j] - substraction;
            }
        }

        Matrix vectorsZeroMean = new Matrix(result);

        //Calculate the covariance matrix V
        //Calculate the square of zeroMean
        double[][] res = new double[m][n];
        for (int i = 0; i < m; ++i) {
            res[i][i] = 0.0;
            for (int k = 0; k < n; ++k) {
                res[i][i] += result[i][k] * result[i][k];
            }
            for (int j = 0; j < i; ++j) {
                res[i][j] = 0.0;
                for (int k = 0; k < n; ++k) {
                    res[i][j] += result[i][k] * result[j][k];
                }
                res[j][i] = res[i][j];
            }
        }

        Matrix squareZeroMean = new Matrix(res);
        covarianceMatrix = square(vectorsZeroMean).times(1.0 / n);

        //Calculate the eigenvalues and eigenvectors.
        EigenvalueDecomposition decomposition = new EigenvalueDecomposition(covarianceMatrix);
        Matrix dEigenValues = decomposition.getD();
        double[] eigenValues = new double[dEigenValues.getRowDimension()];
        for (int i = 0; i < dEigenValues.getRowDimension(); i++) {
            eigenValues[i] = dEigenValues.get(i, i);
        }
        Matrix eigenVectors = decomposition.getV();

        //Calculate the whitening matrix.
        double[] tEigenValues = invVector(sqrtVector(eigenValues));
        Matrix diagonalEv = new Matrix(tEigenValues.length, tEigenValues.length);
        for (int i = 0; i < tEigenValues.length; i++) {
            diagonalEv.set(i, i, tEigenValues[i]);
        }

        //Matrix whiteningMatrix = diagonalEv.times(eigenVectors.transpose());
        Matrix diag = new Matrix(diag(invVector(sqrtVector(eigenValues))));
        whiteningMatrix = diag.times(eigenVectors.transpose());

        //Finally, calculate the whitened data with the whitening matrix.
        Matrix whitenedData = whiteningMatrix.times(vectorsZeroMean);

        return whitenedData.getArrayCopy();
    }

    /**
     * Inverts every element of the vector.
     * @param inVector the vector
     * @return the resulting vctor
     */
    private static double[] invVector(
            double[] inVector) {
        int m = inVector.length;
        double[] outVector = new double[m];
        for (int i = 0; i < m; ++i) {
            outVector[i] = 1 / inVector[i];
        }
        return (outVector);
    }

    /**
     * 
     * @param diag
     * @return 
     */
    public static double[][] diag(double[] diag) {
        int n = diag.length;
        double[][] res = new double[n][n];
        for (int i = 0; i < n; ++i) {
            res[i][i] = diag[i];
        }
        return (res);
    }

    /**
     * Square roots every element of the vector.
     * @param inVector the vector
     * @return the resulting vector
     */
    private static double[] sqrtVector(
            double[] inVector) {
        int m = inVector.length;
        double[] outVector = new double[m];
        for (int i = 0; i < m; ++i) {
            outVector[i] = Math.sqrt(inVector[i]);
        }
        return (outVector);
    }

    /**
     * Calculates the square matrix A * A'.
     * @param mat
     * @return 
     */
    public static Matrix square(
            Matrix matrix) {
        double[][] mat = matrix.getArrayCopy();
        int m = mat.length;
        int n = mat[0].length;
        double[][] res = new double[m][m];
        for (int i = 0; i < m; ++i) {
            res[i][i] = 0.0;
            for (int k = 0; k < n; ++k) {
                res[i][i] += mat[i][k] * mat[i][k];
            }
            for (int j = 0; j < i; ++j) {
                res[i][j] = 0.0;
                for (int k = 0; k < n; ++k) {
                    res[i][j] += mat[i][k] * mat[j][k];
                }
                res[j][i] = res[i][j];
            }
        }
        return new Matrix(res);
    }

    /**
     * 
     * @param inMatrix
     * @param power
     * @return
     */
    private static Matrix powerSymmMatrix(
            Matrix inMatrix,
            double power) {
        EigenvalueDecomposition eigenDeco =
                new EigenvalueDecomposition(inMatrix);
        int m = inMatrix.getRowDimension();
        double[][] eigenVectors = eigenDeco.getV().getArray();
        double[][] eigenValues = eigenDeco.getD().getArray();
        for (int i = 0; i < m; ++i) {
            eigenValues[i][i] = Math.pow(eigenValues[i][i], power);
        }
        Matrix eVectors = new Matrix(eigenVectors);
        Matrix eValues = new Matrix(eigenValues);
        return eVectors.times(eValues).times(eVectors.transpose());
    }

    /**
     * Performs Blind Source Seperation on an ECG Signal with multiple channels.
     * @param source
     * @return 
     */
    public double[][] BSS(double[][] source, int maxIterations) {
        int channels = source.length, m = source.length, noComponents = source.length;
        int samples = source[0].length, n = source[0].length;
        Matrix sourceMatrix = new Matrix(source);
        double[][] component = new double[noComponents][samples];

        //Setup the bBar Matrix
        Matrix bBarMatrix = this.random;
        bBarMatrix = powerSymmMatrix(square(bBarMatrix), -0.5).times(bBarMatrix);
        double epsilon = 0.001;

        for(int k = 0; k < maxIterations; k++) {
            Matrix oldbBar = bBarMatrix.copy();
            for (int curComponent = 0; curComponent < noComponents; curComponent++) {
                //Step 1
                //Create initial vector w.
                Matrix w = oldbBar.getMatrix(curComponent, curComponent, 0, bBarMatrix.getColumnDimension()-1);

                //Step 2
                //Expectation part.
                Matrix sum = new Matrix(channels, 1);
                for (int i = 0; i < samples; i++) {
                    Matrix res1 = w.times(sourceMatrix.getMatrix(0, channels - 1, i, i));
                    double scalar = res1.get(0, 0);
                    scalar = Math.pow(scalar, 3);
                    Matrix result = sourceMatrix.getMatrix(0, channels - 1, i, i).times(scalar);
                    sum = sum.plus(result);
                }

                double nElements = 1.0 / samples;
                Matrix expectation = sum.times(nElements);

                //Second part of calculation.
                Matrix w3 = w.transpose().times(3.0);
                Matrix newW = expectation.minus(w3);

                bBarMatrix.setMatrix(curComponent, curComponent, 0, bBarMatrix.getColumnDimension() - 1, newW.transpose());
            }
            //Step 3
            bBarMatrix = powerSymmMatrix(square(bBarMatrix), -0.5).times(bBarMatrix);

            //Step 4
            double matrixDelta = deltaMatrices(bBarMatrix, oldbBar);
            //System.out.println("Matrix delta: " + matrixDelta);
            if (matrixDelta < epsilon) {
                //System.out.println("Converged after " + k + " iterations.");
                break;
            }


        }

        Matrix result = bBarMatrix.times(whiteningMatrix);
        return result.times(sourceMatrix).getArray();
    }

    /**
     * Calculates a difference measure of two matrices
     * relative to their size.
     * @param mat1 the first matrix
     * @param mat2 the second matrix
     * @return the difference measure
     */
    private static double deltaMatrices(
            Matrix mat1,
            Matrix mat2) {
        Matrix test = mat1.minus(mat2);
        double delta = 0.0;
        int m = mat1.getRowDimension();
        int n = mat2.getColumnDimension();
        for (int i = 0; i < m; ++i) {
            for (int j = 0; j < n; ++j) {
                delta += Math.abs(test.get(i, j));
            }
        }
        return (delta / (m * n));
    }
}
